#include "grav.h"
#include "background.h"

extern output op;  // global variable

void ini_grav_params( grav_params *gp, dictionary *dict ){
    ini_cosmo_params( &(gp->cp), dict );

    gp->a_ini           = iniparser_getdouble(dict,"grav:a_ini", 1e-8);
    gp->h_output_size   = iniparser_getint(dict,"grav:h_output_size", 10000);

/*
    printf("gp->a_ini = %g\n", gp->a_ini);
    printf("gp->h_output_size = %d\n", gp->h_output_size);
    exit(0);
*/
}

/**********************************************************************
    y[0] = h
    y[1] = h'
    y[2] = a
    
    ODEs:
    dy[0]/deta = y[1]
    dy[1]/deta = -2*y[2]*H(a)*y[1] - k*k*y[0]
    dy[2]/deta = y[2]*y[2]*H(a)/c_light

    Notice: eta has unit of Mpc
***********************************************************************/

int ODE_grav( double eta, const double y[], double f[], void *params ){

    grav_params *gp = (grav_params*)params;

    double h    = y[0];
    double h1   = y[1];
    double a    = y[2];
    double k    = gp->k;    // Mpc^{-1}

    double Ha   = H_a( a, &(gp->cp) ); // km/s/Mpc
    
    double c_light = 3e5;   // km/s
    
    f[0] = h1;
    f[1] =-2.0*a*Ha*h1/c_light - k*k*h;
    f[2] = a*a*Ha/c_light;
    
    return GSL_SUCCESS;
}

void do_grav_calculation( grav_params *gp, double *k_array, int k_num ){

    int ndim = 3;
    gsl_odeiv2_system   sys = { ODE_grav, NULL, ndim, gp};
    gsl_odeiv2_driver   * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd, 1e-6, 1e-6, 0.0);

    char fout[1024];

    /* Loop over k-modes */
    
    int k_cnt = 0;
    for( k_cnt = 0; k_cnt < k_num; k_cnt++ ){

        double k = k_array[k_cnt];   // for test right now ...
        
        char fname[128];
        sprintf(fname, "k=%d.txt", k_cnt);

        strcpy(fout, op.out_dir);
        strcat(fout, "/");
        strcat(fout, fname);
        
        printf("fout --> %s\n", fout);
        
        FILE *fp = fopen(fout, "w");
        
        gp->k = k;

        /*************************************************************/
        /*  setup initial condition */
        double c_light  = 3e5;   // km/s
        double a_ini    = gp->a_ini;
        double eta_ini  = a_ini*c_light;    // Why?
        double h_ini    = 1.0;
        double h1_ini   =-0.5*eta_ini*k*k*h_ini;
        /*************************************************************/
        
        double eta  = eta_ini;
        double y[3] = {h_ini, h1_ini, a_ini};

        double eta_max = 14551; // this precalculated ...
        double lneta_ini = log(eta_ini);
        double lneta_max = log(eta_max);
        double dlneta = (lneta_max-lneta_ini)/(gp->h_output_size-1.);
        
        fprintf(fp, "#  eta   h    h'    a\n");
        fprintf(fp,"%15.12f %15.12f %15.12f %15.12f\n", eta, y[0], y[1], y[2]);
        
        int i;
        for( i=1; i<gp->h_output_size; i++){
            
            double eta_next = exp( lneta_ini + i*dlneta );

            int status = gsl_odeiv2_driver_apply (d, &eta, eta_next, y);

            if( status != GSL_SUCCESS ){
                printf("failed to evolve ODE ...\n");
                break;
            }

            fprintf(fp,"%15.12f %15.12f %15.12f %15.12f\n", eta_next, y[0], y[1], y[2]);
        }
        fclose(fp);
    }

    gsl_odeiv2_driver_free(d);
        
}
